library(survey)
library(foreign)
library(ggplot2)
library(ggrepel)
library(cowplot)
library(dplyr)
library(glmnet)

## see PCM_cleaning.R for data cleaning
#load("~/Dropbox/Public_Conf_Mil/Data_and_code/wave2_data/clean_dataw2.RData")
# This is saved: C:\Users\Robert Allred\Documents\Duke\Civ_Mil_Appendix\Code

## Create weighted survey design object
w2_design <-
  svydesign(
    id = ~ 1,
    weights = ~ weight2,
    data = df
  )

#### FIGURE 2.9 - DEMOGRAPHIC PREDICTIONS ####

demo_mod <- svyglm(Q11_d ~ dem + rep + ideo3 + male + white + black + hispanic + asian + income5 +
                     EDUC5 + unemployed + silent + boomer + genx + millen + activeduty + vet +
                     social + family + midwest + south + west + catholic + christian + norelig +
                     city + rural + married + A2 + A3 + A4 + A5 + A6 + A7 + A8,
                   family = binomial(link = "logit"),
                   design = w2_design)
summary(demo_mod)




### Predictions by Demographic Profiles

## Demo 1: Asian female student
demo1 <- data.frame(dem = 1, rep = 0, ideo3 = 0,
                    male = 0, white = 0, black = 0, hispanic = 0, asian = 1,
                    income5 = 1, EDUC5 = 3, unemployed = 0,
                    silent = 0, boomer = 0, genx = 0, millen = 0,
                    activeduty = 0, vet = 0, social = 0, family = 0,
                    midwest = 0, south = 0, west = 1,
                    catholic = 0, christian = 0, norelig = 1,
                    city = 0, rural = 0, married = 0,
                    A2 = 0, A3 = 0, A4 = 0, A5 = 0, A6 = 0, A7 = 0, A8 = 0)
demo1pred <- as.data.frame(predict(demo_mod, newdata = demo1, type = "response", se.fit = TRUE))

## Demo 2: Black low-income single female
demo2 <- data.frame(dem = 1, rep = 0, ideo3 = 1,
                    male = 0, white = 0, black = 1, hispanic = 0, asian = 0,
                    income5 = 2, EDUC5 = 2, unemployed = 0,
                    silent = 0, boomer = 0, genx = 0, millen = 1,
                    activeduty = 0, vet = 0, social = 0, family = 0,
                    midwest = 0, south = 1, west = 0,
                    catholic = 0, christian = 0, norelig = 1,
                    city = 1, rural = 0, married = 0,
                    A2 = 0, A3 = 0, A4 = 0, A5 = 0, A6 = 0, A7 = 0, A8 = 0)
demo2pred <- as.data.frame(predict(demo_mod, newdata = demo2, type = "response", se.fit = TRUE))

## Demo 3: CA white collar male worker
demo3 <- data.frame(dem = 0, rep = 0, ideo3 = 1,
                    male = 1, white = 1, black = 0, hispanic = 0, asian = 0,
                    income5 = 5, EDUC5 = 4, unemployed = 0,
                    silent = 0, boomer = 0, genx = 0, millen = 1,
                    activeduty = 0, vet = 0, social = 0, family = 0,
                    midwest = 0, south = 0, west = 1,
                    catholic = 0, christian = 0, norelig = 1,
                    city = 1, rural = 0, married = 0,
                    A2 = 0, A3 = 0, A4 = 0, A5 = 0, A6 = 0, A7 = 0, A8 = 0)
demo3pred <- as.data.frame(predict(demo_mod, newdata = demo3, type = "response", se.fit = TRUE))

## Demo 4: Hispanic male migrant worker
demo4 <- data.frame(dem = 0, rep = 0, ideo3 = 1,
                    male = 1, white = 0, black = 0, hispanic = 1, asian = 0,
                    income5 = 1, EDUC5 = 1, unemployed = 0,
                    silent = 0, boomer = 0, genx = 1, millen = 0,
                    activeduty = 0, vet = 0, social = 0, family = 0,
                    midwest = 0, south = 0, west = 1,
                    catholic = 1, christian = 0, norelig = 0,
                    city = 0, rural = 1, married = 1,
                    A2 = 0, A3 = 0, A4 = 0, A5 = 0, A6 = 0, A7 = 0, A8 = 0)
demo4pred <- as.data.frame(predict(demo_mod, newdata = demo4, type = "response", se.fit = TRUE))

## Demo 5: CA liberal male PhD
demo5 <- data.frame(dem = 1, rep = 0, ideo3 = 0,
                    male = 1, white = 1, black = 0, hispanic = 0, asian = 0,
                    income5 = 4, EDUC5 = 5, unemployed = 0,
                    silent = 0, boomer = 0, genx = 1, millen = 0,
                    activeduty = 0, vet = 0, social = 0, family = 0,
                    midwest = 0, south = 0, west = 1,
                    catholic = 0, christian = 0, norelig = 1,
                    city = 1, rural = 0, married = 0,
                    A2 = 0, A3 = 0, A4 = 0, A5 = 0, A6 = 0, A7 = 0, A8 = 0)
demo5pred <- as.data.frame(predict(demo_mod, newdata = demo5, type = "response", se.fit = TRUE))

## Demo 6: NE Boomer woman
demo6 <- data.frame(dem = 1, rep = 0, ideo3 = 1,
                    male = 0, white = 1, black = 0, hispanic = 0, asian = 0,
                    income5 = 4, EDUC5 = 3, unemployed = 0,
                    silent = 0, boomer = 1, genx = 0, millen = 0,
                    activeduty = 0, vet = 0, social = 0, family = 0,
                    midwest = 0, south = 0, west = 0,
                    catholic = 0, christian = 0, norelig = 1,
                    city = 0, rural = 0, married = 1,
                    A2 = 0, A3 = 0, A4 = 0, A5 = 0, A6 = 0, A7 = 0, A8 = 0)
demo6pred <- as.data.frame(predict(demo_mod, newdata = demo6, type = "response", se.fit = TRUE))

## Demo 7: Midwestern soccer Mom
demo7 <- data.frame(dem = 0, rep = 0, ideo3 = 1,
                    male = 0, white = 1, black = 0, hispanic = 0, asian = 0,
                    income5 = 4, EDUC5 = 3, unemployed = 0,
                    silent = 0, boomer = 0, genx = 1, millen = 0,
                    activeduty = 0, vet = 0, social = 0, family = 1,
                    midwest = 1, south = 0, west = 0,
                    catholic = 0, christian = 1, norelig = 0,
                    city = 0, rural = 0, married = 1,
                    A2 = 0, A3 = 0, A4 = 0, A5 = 0, A6 = 0, A7 = 0, A8 = 0)
demo7pred <- as.data.frame(predict(demo_mod, newdata = demo7, type = "response", se.fit = TRUE))

## Demo 8: Black female NCO
demo8 <- data.frame(dem = 1, rep = 0, ideo3 = 0,
                    male = 0, white = 0, black = 1, hispanic = 0, asian = 0,
                    income5 = 3, EDUC5 = 2, unemployed = 0,
                    silent = 0, boomer = 0, genx = 0, millen = 1,
                    activeduty = 1, vet = 0, social = 1, family = 1,
                    midwest = 0, south = 1, west = 0,
                    catholic = 0, christian = 1, norelig = 0,
                    city = 1, rural = 0, married = 0,
                    A2 = 0, A3 = 0, A4 = 0, A5 = 0, A6 = 0, A7 = 0, A8 = 0)
demo8pred <- as.data.frame(predict(demo_mod, newdata = demo8, type = "response", se.fit = TRUE))

## Demo 9: Post-911 white male veteran
demo9 <- data.frame(dem = 0, rep = 1, ideo3 = 2,
                    male = 1, white = 1, black = 0, hispanic = 0, asian = 0,
                    income5 = 3, EDUC5 = 3, unemployed = 0,
                    silent = 0, boomer = 0, genx = 0, millen = 1,
                    activeduty = 0, vet = 1, social = 1, family = 0,
                    midwest = 0, south = 0, west = 1,
                    catholic = 1, christian = 0, norelig = 0,
                    city = 0, rural = 1, married = 0,
                    A2 = 0, A3 = 0, A4 = 0, A5 = 0, A6 = 0, A7 = 0, A8 = 0)
demo9pred <- as.data.frame(predict(demo_mod, newdata = demo9, type = "response", se.fit = TRUE))

## Demo 10: Hispanic male marine
demo10 <- data.frame(dem = 0, rep = 0, ideo3 = 0,
                     male = 1, white = 0, black = 0, hispanic = 1, asian = 0,
                     income5 = 3, EDUC5 = 3, unemployed = 0,
                     silent = 0, boomer = 0, genx = 1, millen = 0,
                     activeduty = 1, vet = 0, social = 1, family = 0,
                     midwest = 0, south = 0, west = 1,
                     catholic = 1, christian = 0, norelig = 0,
                     city = 1, rural = 0, married = 0,
                     A2 = 0, A3 = 0, A4 = 0, A5 = 0, A6 = 0, A7 = 0, A8 = 0)
demo10pred <- as.data.frame(predict(demo_mod, newdata = demo10, type = "response", se.fit = TRUE))

## Demo 11: NASCAR Dad
demo11 <- data.frame(dem = 0, rep = 1, ideo3 = 2,
                     male = 1, white = 1, black = 0, hispanic = 0, asian = 0,
                     income5 = 3, EDUC5 = 3, unemployed = 0,
                     silent = 0, boomer = 0, genx = 1, millen = 0,
                     activeduty = 0, vet = 0, social = 1, family = 1,
                     midwest = 0, south = 1, west = 0,
                     catholic = 0, christian = 1, norelig = 0,
                     city = 0, rural = 1, married = 1,
                     A2 = 0, A3 = 0, A4 = 0, A5 = 0, A6 = 0, A7 = 0, A8 = 0)
demo11pred <- as.data.frame(predict(demo_mod, newdata = demo11, type = "response", se.fit = TRUE))

## Demo 12: NE white collar worker
demo12 <- data.frame(dem = 0, rep = 1, ideo3 = 1,
                     male = 1, white = 1, black = 0, hispanic = 0, asian = 0,
                     income5 = 4, EDUC5 = 4, unemployed = 0,
                     silent = 0, boomer = 0, genx = 1, millen = 0,
                     activeduty = 0, vet = 0, social = 0, family = 0,
                     midwest = 0, south = 0, west = 0,
                     catholic = 0, christian = 1, norelig = 0,
                     city = 1, rural = 0, married = 1,
                     A2 = 0, A3 = 0, A4 = 0, A5 = 0, A6 = 0, A7 = 0, A8 = 0)
demo12pred <- as.data.frame(predict(demo_mod, newdata = demo12, type = "response", se.fit = TRUE))

## Demo 13: Male Midwestern farmer
demo13 <- data.frame(dem = 0, rep = 1, ideo3 = 2,
                     male = 1, white = 1, black = 0, hispanic = 0, asian = 0,
                     income5 = 2, EDUC5 = 2, unemployed = 0,
                     silent = 0, boomer = 1, genx = 0, millen = 0,
                     activeduty = 0, vet = 0, social = 1, family = 0,
                     midwest = 1, south = 0, west = 0,
                     catholic = 1, christian = 0, norelig = 0,
                     city = 0, rural = 1, married = 1,
                     A2 = 0, A3 = 0, A4 = 0, A5 = 0, A6 = 0, A7 = 0, A8 = 0)
demo13pred <- as.data.frame(predict(demo_mod, newdata = demo13, type = "response", se.fit = TRUE))

## Demo 14: Female Fox News watcher
demo14 <- data.frame(dem = 0, rep = 1, ideo3 = 2,
                     male = 0, white = 1, black = 0, hispanic = 0, asian = 0,
                     income5 = 3, EDUC5 = 2, unemployed = 0,
                     silent = 0, boomer = 1, genx = 0, millen = 0,
                     activeduty = 0, vet = 0, social = 1, family = 1,
                     midwest = 0, south = 0, west = 0,
                     catholic = 0, christian = 1, norelig = 0,
                     city = 0, rural = 0, married = 1,
                     A2 = 0, A3 = 0, A4 = 0, A5 = 0, A6 = 0, A7 = 0, A8 = 0)
demo14pred <- as.data.frame(predict(demo_mod, newdata = demo14, type = "response", se.fit = TRUE))

## Demo 15: Male Republican Korean war veteran
demo15 <- data.frame(dem = 0, rep = 1, ideo3 = 2,
                     male = 1, white = 1, black = 0, hispanic = 0, asian = 0,
                     income5 = 4, EDUC5 = 3, unemployed = 0,
                     silent = 0, boomer = 1, genx = 0, millen = 0,
                     activeduty = 0, vet = 1, social = 1, family = 1,
                     midwest = 0, south = 0, west = 1,
                     catholic = 0, christian = 1, norelig = 0,
                     city = 0, rural = 0, married = 0,
                     A2 = 0, A3 = 0, A4 = 0, A5 = 0, A6 = 0, A7 = 0, A8 = 0)
demo15pred <- as.data.frame(predict(demo_mod, newdata = demo15, type = "response", se.fit = TRUE))

demo_preds <- rbind(demo2pred, demo1pred, demo5pred, demo3pred, 
                    demo6pred, demo7pred, demo4pred, demo12pred, demo8pred, 
                    demo10pred, demo14pred, demo13pred, demo11pred, 
                    demo9pred, demo15pred)
demo_preds$ci_lo <- demo_preds$response - 1.96*demo_preds$SE
demo_preds$ci_hi <- demo_preds$response + 1.96*demo_preds$SE
demo_preds$count <- c(1:15)

demo_plot <- ggplot(demo_preds) +
  geom_pointrange(aes(x = count, y = response, ymin = ci_lo, ymax = ci_hi)) +
  coord_flip() + xlab("") + ylab("Predicted Support for the Military") +
  scale_x_continuous(breaks = c(1:15),
                     labels = c("Black low-income single female", 
                                "Asian female student", 
                                "CA liberal male PhD", 
                                "CA male white collar worker",
                                "NE Boomer woman", 
                                "Midwestern soccer mom", 
                                "Hispanic male migrant worker", 
                                "NE white collar worker",
                                "Black female NCO", 
                                "Hispanic male marine",
                                "Female Fox News watcher",
                                "Male Midwestern farmer",
                                "NASCAR Dad",
                                "Post-9/11 white male veteran",
                                "Male Republican Korean war veteran")) +
  theme_bw()
demo_plot
# ggsave("~/Dropbox/Public_Conf_Mil/Chapter 2 Figures/fig2-10_final.jpeg", plot = demo_plot, dpi = 1000)
